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Abstract. The first passage time density of a diffusion process to a time 
varying threshold is of primary interest in different fields. Here, we consider a 
Brownian motion in presence of an exponentially decaying threshold to model 
the neuronal spiking activity. Since analytical expressions of the first passage 
time density are not available, we propose to approximate the curved boundary 
by means of a continuous two-piecewise linear threshold. Explicit expressions 
for the first passage time density towards the new boundary are provided. First, 
we introduce different approximating linear thresholds. Then, we describe how 
to choose the optimal one minimizing the distance to the curved boundary, and 
hence the error in the corresponding passage time density. Theoretical means, 
variances and coefficients of variation given by our method are compared with 
empirical quantities from simulated data. Moreover, a further comparison with 
firing statistics derived under the assumption of a small amplitude of the time- 
dependent change in the threshold, is also carried out. Finally, maximum 
likelihood and moment estimators of the parameters of the model are derived 
and applied on simulated data. 


1. Introduction. Stochastic models have been extensively used in theoretical neu¬ 
roscience since the pioneer work by Gerstein and Mandelbrot in 1964 [12], There 
they considered a Wiener process (also known as Brownian motion or Perfect- 
Integrate-and-Fire model) to model the voltage across the membrane. An action po¬ 
tential, also known as spike, is generated whenever the membrane potential reaches 
a certain constant threshold. After that, the membrane voltage is reset to its rest¬ 
ing value and the evolution restarts. From a mathematical point of view, a spike is 
the first passage time (FPT) of a stochastic process to a constant threshold. The 
collection of spike epochs of a neuron, called spike train, defines a renewal process, 
with independent and identically distributed inter-spike intervals (ISIs). Despite 
the excellent fit with some experimental data, the Gerstein-Mandelbrot model was 
criticized because it disregards features involved in neuronal coding. 
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A first extension, combining both mathematical tractability and biological re¬ 
alism, is represented by Leaky-Integrate-and-Fire (LIF) models [28, 36]. Despite 
some criticisms on the lack of fit of experimental data [16, 31], these models are still 
largely used. 

Another common generalization is represented by Wiener processes (or more 
generally LIF models) with time-dependent threshold [35, 37]. These models can be 
chosen to reproduce biological features such as the afterhyperpolarization in neu¬ 
rons. For exponentially decaying thresholds, these processes can be used to model a 
neuron with an exponential time-dependent drift, as shown by Lindner and Longtin 
[19]. They investigated the effect of an exponentially decaying threshold on the br¬ 
ing statistics of a stochastic integrate-and-fire neuron [19]. Using a perturbation 
method [18], they derived analytical expressions of the firing statistics under the 
assumption that the amplitude e of the time-dependent change in the threshold is 
small. These statistics are useful to characterize the spontaneous neural activity 
and to investigate the neuronal signal transmission. In particular, they can suggest 
under which conditions a decaying threshold may facilitate or deteriorate signal 
processing by stochastic neurons. For a Wiener process, these quantities can also 
be obtained using the approach in [38]. Also this method assumes a small ampli¬ 
tude e, but it has the advantage of providing an explicit approximation of the FPT 
density. 

Here we consider a Wiener process with exponentially decaying threshold. The 
first aim of the paper is to provide an alternative method to approximate the firing 
statistics and the FPT density for any possible amplitude e, extending the results 
in [19, 38]. Different estimators are proposed, as mentioned in Section 1.1 and 
discussed in Section 4.2. Means, variances, coefficients of variation (CVs) and dis¬ 
tributions of the FPTs are compared on simulated data and the most suitable are 
recommended. A comparison with the results in [19, 38] under the assumption of a 
small amplitude e is also performed. The second aim of this work is the estimation 
of drift and diffusion coefficients of the Wiener process. Maximum likelihood and 
moment estimators are derived and evaluated on simulated data. Our results show 
a good approximation of both firing statistics and parameters of the underlying 
model. 

Although the considered model generates a renewal process, the proposed method 
can also be applied to non-renewal processes, e.g. adaptive threshold models [8, 17]. 
Recently, an increasing interest arose towards these models, interest motivated by 
the excellent fit of the firing statistics of electrosensory neurons [7, 9] . The novelty 
of these models is that the threshold has a jump immediately after a spike. Since 
the boundary depends on the previous firing epochs, the ISIs are not independent 
anymore. However, the distribution between two consecutive spikes, conditioned on 
the initial position of the threshold, is the same of that studied here. Hence, our 
results may represent a first step towards an understanding of the more complicated 
adapting-threshold models. 

1.1. Mathematical background. FPTs of diffusion processes to constant or time- 
dependent thresholds have been extensively studied in the literature. Explicit ex¬ 
pressions for constant thresholds are available for the Wiener process [10, 11], for 
a special case of the Ornstein Uhlenbeck (OU) process [26], for the Cox-Ingersoll- 
Ross process [6], and for those processes which can be obtained from the previous 
through suitable measure or space-time transformations, see e.g. [2, 6, 25]. For most 
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of the processes arising from applications and for time-varying thresholds, analyt¬ 
ical expressions are not available. Numerical algorithms based on solving integral 
equations have been proposed in [4, 5, 27, 29, 32, 34], while approximations based 
on Monte-Carlo path-simulation methods in [13, 14, 20]. 

A different approach to tackle the FPT problem consists in focusing directly 
on the two-sided boundary crossing probability (BCP), i.e. the probability that 
a process is constrained to be between two boundaries. If one of the boundary 
is set to —oo, the resulting one-sided BCP equals the survival probability of the 
FPT to the other boundary [39]. Explicit formulas for the BCP of a standard 
Wiener process for continuous and piecewise-linear thresholds are known (see [3, 
22, 23, 39, 40]). In general, the BCP of a diffusion process through an exponential 
decaying threshold is available only for those processes which can be expressed as 
a piecewise monotone functional of a standard Brownian motion. Examples are 
the OU process or the geometric Brownian motion with time dependent drift for 
specific parameter values [40] . The simple but powerful idea is to approximate both 
one and two-sided curvilinear BCPs by similar probabilities for close boundaries of 
simpler form, namely n piecewise-linear thresholds, whose computation of the BCP 
for Wiener is feasible. Under some mild assumptions, the approximated two-sided 
BCP converges to the original one when n —^ oo [40] , with rate of convergence given 
in [3, 39]. 

For the exponential decaying threshold considered in this paper, the conver¬ 
gence can be obtained by choosing piecewise linear thresholds approximating the 
curved boundary from above and below, with approximation accuracy given by 
their distance [40]. However, all the available formulas for the BCPs require either 
Monte-Carlo simulation methods or heavy numerical approximations. 

Here we consider a two-piecewise linear threshold as an approximation of the 
curvilinear boundary. Since n = 2, the asymptotic convergence of the BCPs does 
not hold. However, we can derive analytical expression for the FPT density to the 
two-piecewise linear boundary, and use it as an approximation of the unknown FPT 
density. Four possible piecewise thresholds are proposed and optimized to minimize 
the distance to the original threshold. 

2. Model. We describe the membrane potential evolution of a single neuron by a 
Wiener process X{t), starting at some initial value xq- We assume X{t) given as 
the solution to a stochastic differential equation 

f dX(t) = ^dt + adW(t), , , 

\ X{to) =xo, t> to, 

where W{t) is a standard (driftless) Wiener process. The drift /r > 0 and the 
diffusion coefficient ct > 0 represent input and noise intensities, respectively. A 
spike occurs when the membrane potential X (t) exceeds the exponentially decaying 
threshold 

6*(t) = 6o-I-eexp [-A(t - 4)] (2) 

for the first time. Here, Sk denotes the time of the fcth spike for fc > 0, and can be 
interpreted as a relative refractory period. We set 4 to be the starting time of the 
process, i.e. 4 = 4- The term A represents the decay rate of the threshold, while 
e is interpreted as the amplitude of the time-dependent change in the boundary. 
After a spike, the membrane potential is reset to its resting position xq < bo + e, 
and its evolution is restarted, as illustrated in Fig. I. The presence of 4 in (2) 
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ensures that the ISIs are independent and identically distributed. Denote by b{t) 
the threshold b*{t) for fc = 0, i.e. 

b{t) =bo + e exp [-\{t - to)] • 

Then, all ISIs are distributed as the FPT of X to b{t), namely 

Tb = inf{t > to : X{t) > b{t)}. 

Quantities of interest are the probability density function (pdf) and the cumulative 
distribution function (cdf) of Tt, denoted by and Ft^, respectively. Another 
relevant quantity is the two-sided BCP given by 

Px(a, c, r) = P(o(t) < X{t) < c(t),Vt S [to,T]). 

Here t > to is fixed, boundaries a{t) and c(t) are real functions satisfying a{t) < c(t) 
for Si\\to<t<T and a{to) < Xo < c{to). Setting a{t) = —oo yields the one-sided 
BCP 

Px(-oo,c,r) = P(Tc > r) = 1 - 

which corresponds to the survival probability of Tc- For a standard Wiener process 
W, Wang and Potzelberger [40] showed that, if the sequences of piecewise linear 
functions a„ and c„ converge uniformly to a{t) and c(t) on [toi^] respectively, then, 
for the continuity property of probability measure, it holds 

lim Pwian,Cn,T) = Pw{a,c,T). 

n—>-oo 

When a{t) = —oo and c{t) = b{t), the convergence of P(rb„ > t) to P(Tb > r) can 
be obtained by choosing piecewise linear thresholds approximating b{t) from above, 
denoted by b^{t), or from below, b~(t). That is, 6+(t) > b{t) and b~{t) < b{t), Vt S 
[to,T], respectively. Since the considered curved boundary is convex, we have [23] 

f‘x{-co,b~,T) < Px(-oo,6,t) < P(-oo,6+,r), (3) 

i.e. 

P(T^+ < r) < P(rb < r) < P(r^- < t). 

The approximation accuracy is given by Pv:(—oo, t)—P x(—oo, b~,T) = Ft^_ (j)— 
with bounds given in [3]. Obviously, the accuracy in the BCP increases 
when the distance between the two thresholds decreases. 


3. FPT to continuous piecewise linear threshold. The transition density 
function of a standard Brownian motion in xi, X 2 ,..., x„ at time ti, O, • ■ •, tn, con¬ 
strained to be below the absorbing threshold c{t) defined by n piecewise-linear 
threshold over [to,tn]j is given in [39]. Extending that result to the case of a Brow¬ 
nian motion with drift /r and diffusion coefficient cr, starting in xq < c(to) = Cq at 
time toj we obtain 


n 

Pc • ■ • : Pc 5 | 17 1 ) 

2=1 


= n 


[1 - exp 




A/27rcr2(t- - ti_i) 


exp 


[Xi - Xi i 


20-2 (i- 


- U-i) 


,(4) 


where ci = c(ti) and Xi < Ci,l < i < n. 
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Figure 1 . Schematic illustration of the single trial of a Wiener 
process in presence of exponentially decaying threshold b{t) = bo + 
eexp(—A(t — Sk)), where Sk denotes the A:th spike. The membrane 
potential starts in xq = 0 at time Jq := to = 0, and it evolves until 
it hits the boundary for the first time. Then, a spike is generated, 
the voltage X{t) is reset to its resting potential xq, the threshold 
is reset to 6o + e and the evolution restarts. For the considered 
spiking generation rule, all ISIs are independent and identically 
distributed. Here the parameters are /r = l,cr^ = l,6o = 1;-^ = 1 
and e = 5. 


From (4), it follows that [41] 


¥{Wt, € Cl, ..., € C„, T, > t„) = 


/ -I 

JCi JC, 


Pc{xi,tl- . . .■,X„,tn\xo,to)dXl ■ ■ ■ dXr, 


(5) 


for any Borel set Ci C (— oo,Ci), \ < i < n. If Ci = (—oo,Ci), then (5) is equal to 
P(Tc > tn), and it holds 

d T" 

/r,(t)=-7^/ •■■/ ,Xn,tn\xQ,tQ)dxi---dXn- (6) 

Otji J — oo 


When n = 1, the pdf /t„ is known. Since X is a Wiener process with positive 
drift, the distribution of the FPT to c{t) = a + f3{t — to) is inverse Gaussian, 
Tc ~ IG [(a - xo)/{p - /3), (a - xoY , with pdf 


It, (t) 


a — Xq 

— / exp 

^2'Ka'^{t - to)^ 


Xq - {p - I3){t - to)f 
2a'^{t - to) 


(7) 


mean E[Tc] = (a — xo)/{p — /3) and variance Var(Tc) = {a — xo)(t'^/{ p — /3)^ [10, 
11]. Note that the distribution of is the same of that of the FPT of a Wiener 
process with positive drift ^ — /3 to a constant threshold c{t) = a. In general, the 
approximation of by when n = 1 is too rough. However, when A is very 
small, exp(—At) « 1 —At, yielding 6(t) « 6o + e—Aet. Hence, T}, can be approximated 
by Tc with a = bo + e and (3 = —Ae. 
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Since we approximate b(t) by means of a continuous two-piecewise linear thresh¬ 
old, we denote by b the linear threshold c{t) when n = 2. We have 


b{t) — bi{t)l[t<ti} + &2(t)l{t>ti} — 


o:2 + P2(t-ti) 


( 8 ) 


where 1 a denotes the indicator function of the set A and ai,a 2 ,^i ,/?2 € K- 
Throughout the paper, we set a 2 = ai + /3i(ti — to) to guarantee the continu¬ 
ity of b(t). This allows to provide analytical expressions of (4) and (6), which we 
use as an approximation of /t^ . In particular, we have 

Fm <t) = Fm <t,T~,< h)+Fm <t,T~,> to 

rH^i) 

= < min(ti,t)) -F / P(T^^ < t\X{ti) = xi)pi^{xi,ti)dxi 

J —OO 


= / /tj {s)ds+ / / /tj {s\xi,ti)pi^{xi,ti)dsdxi (9) 

Jo ^ J-oo Jti ^ 

with pg {xi,ti) given by (4) for n = 1. Mimicking [33], by taking the derivative of 
(9) with respect to t, and plugging (7) in it for proper values of a and /3, we get 

rb{ti) 

\ M —/3i ’ J J —OO M —^2 ’ ^ / 


> —OO 

rb(ti) ft 


Ql — Xq 


\/27rcr2(t - toY 


exp 


(oi -Xo- {p- /3i)(t - to))" 




exp 


y/2-na^{t — to)^ 

X Ua2 - xo - P2{ti - to)] 4 > 


20-2 (t - to) 

(q2 - Xq- (p- P2){t - tl) - p{ti - to))^ 

20-2(i - to) 

(02 - 3^0 - P2{tl - to)) V(t - tl) 

o/a^{ti - to){t - to) 

2{t - ti){ai - xo){a 2 - oi - /92(ti - to)) 


( 10 ) 


-(02 + ^0 - Mti - to) - 201) exp 

^ ^ [ (o 2 -f a^o - / 32 (ti - to) - 2 oi)y (t - tl) 

\ - to)(t - to) 

This result extends that for a driftless Brownian motion, see e.g. [1, 30]. As 
expected, setting ai = a 2 = a and Pi = h = jd yields the pdf of the FPT of a 
Wiener process to a linear threshold c(t) = a -\- Pit — to). By definition, the first 
two moments and variance of Tj, are given by 


pOO pOQ 

m]= tfT.m, E[T?]= ffT.m, Var[T^]=E[T|]-E[T^]^ 
Jo Jo 


and can be numerically computed. 


(11) 


4. Parameter estimation. 

4.1. Parameter estimation of the piecewise-linear threshold. The primary 
aim of this paper is the approximation of the FPT distribution (and relevant statis¬ 
tics) for a curved boundary b{t), by means of the FPT distribution for a continuous 
two-piecewise linear threshold bit). As discussed in Section 2, the quality of the 
approximation improves when the distance between b and b decreases. 
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Curved threshold b(t) and two-piecewise threshold b(t) 




Figure 2. Curved threshold b{t) (continuous line) and four pro¬ 
posed approximating piecewise-linear thresholds: 6+ (t) from above 
(dashed lines) ; b-{t) from below (dashed-dotted line); 6betw(i) 
which is equidistant from b^ and b_ (gray dashed line); &free(i) 
with no restrictions (gray continuous line). For each type of linear 
threshold, the best approximation is given by the line minimizing 
a function of the distance from b on [to = 0, r] (left figure) and on 
[to,t*] C [to = 0,r] (right figure). As discussed in Section 4.1, a 
shorter time interval provides a better approximation of b. 


Denote hy 9 = {ai,/3i,/32, ti) the parameters of b in (8), with a 2 = ai-|-/3i(ti—to). 
We are interested in determining the estimator 9 which minimizes \b{t) — b(t)\ on 
[To,r*], with to < To < ti < r* < r. The time interval is chosen such that the 
probability of having a FPT outside it is smaller than 0.005, i.e. 

P(ne [To,r*])>0.99. (12) 

Doing this, we improve the approximation of b on [ro,r*] (cf. Fig. 2), allowing a 
larger deviation from b on [to, To) and (r*,T], i.e. on intervals where the probability 
of observing a FPT is low. Since 

V{Tb >t)= P(A(s) < 6(s),s e [0,t]) < P(A(t) < 6(t)) 

and b{t) > &oi it follows that 

P(A(t) > bit)) < P(Tb < t) < P(T6„ < t), 

with Tbg ~ IG{{bo — xo)/tJ-, {bo — xo)‘^/a^). Since A is a Wiener process, X{t) ^ 
N{fit,a‘^t). Then, we choose tq and r* such that 

P(rf,o < To) = 0.005, P(A(t*) > 5(t,)) = 0.995, 

yielding the desired probability (12). 

Throughout the paper, we consider four possible continuous two-piecewise linear 
boundaries on [to,t*], as illustrated in Fig. 2: 
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1. Threshold 6+ approximating b from above on passing through (tq, 

6(ro)), (ti,6(ti)) and 

b+[t) = b{To) H---- {t - ro)l{t<ii} H- - - {t - 


h - To 


i.e. 


u\ a b{ti)-b{To) ^ 6 (t,) - 

O^l — b+[to), Pi — - 7 -) P2 — - 7 - {t — tl}- 


h - To 


Titi ti 


Due to the assumptions, for given tq and r*, ti is the only unknown quantity. 

2. Threshold approximating b from below on We assume that is 

tangent to b{t) in both ti and t 2 > ti, with ti intersection time point of the 
two tangent lines 

yiit) = b{ii) - Aeexp(-Ati)(t - U), 

for i = l,2. Setting yi{ti) = ?/ 2 (ti), we get 

^ _ exp(-Afi)[l + Ati] - exp(-At 2 )[l + At 2 ] 

A[exp(—Ati) — exp(—At 2 )] 

Then, the desired threshold b-{t) is 

z, zr ^ , Viih) - yiih),, X . J/2(i2)-2/2(ii),, , 

b-{t) = yiih) + ^ X + r ^ 

ti — Cl 12 — Cl 


with 


u u\ p yiih) - yiiii) a ^ 2 (^ 2 )-2/2(ti) 

ai — O_(to), Pi — -;- = -, P2 — 


tl — tl 


t 2 — tl 


For fixed tq and r», the unknown parameters are ti and ^ 2 - 

3. Threshold &betw(t) constrained to be between 6+(t) and b-{t) on [ro,r*], i.e. 
b-{t) < bbetwit) < b+it). 

4. Threshold with no constraints, denoted by bfree{t). 

Denote by 6*+, 0_, 0betw and 6*free the estimators of 0 from the boundaries 6+, 6_, 6betw 
and &free, respectively. From (3), it follows that the best approximation of Pv:(— 00 , b, 
t) is obtained when the distance between &_|_ and 6_ is minimized. For this reason, 
we define 0+ and 0_ as the estimators minimizing the area of the squared distance 
between the two boundaries, i.e. 


(0+, ) = arg min 


\b+{t) — b-{t)\^dt 


with 0+ and 9_ satisfying the conditions 6+(t) > b{t) and 6_(t) < b{t) on [ro,r*]. 
The quantity |6+(t) — b_{t)\‘^ instead of |6+(t) — b_{t)\ is chosen to avoid possible 
numerical issues in the optimization procedure. 

Once 6+ and 6_ have been computed, the estimator 0betw is defined as 


6 'betw = arg min 


{\b+{t) - 6 betw(t)P + \b-{t) - 6 betw(t)P) dt 


i.e. it is the equidistant line from 6+(t) and b-(t). 
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Finally, the estimator is the one minimizing the area of the squared distance 
between the piecewise-line and the curved boundary, i.e. 

[ \bireeit) -b(t)\'^dt . 

To 

Note that the estimation of 6 does not depend on the observations of the FPTs but 
it is performed theoretically under the assumption that the parameters A, e and bg of 
b{t) are known. The proposed estimators and their assumptions are summarized in 
Table 1. All the minimizations have been performed in the computing environment 
R [24]. Since the parameter values need to fulfil some conditions (cf. Table 1), 
minimizing the areas is a constrained optimization problem. We perform it by 
means of the built-in R function optim, penalizing those parameter values not 
fulfilling the conditions by returning 10^°. 

4.2. Parameter estimation of the process. The second aim of the paper is the 
estimation of the parameters p, and cr^ of the Wiener process from a sample 
of n independent observations of Tf,. That is, we want to estimate (j) = {p, a^) under 
the assumption that the parameters of the threshold are known. 

4.2.1. Maximum likelihood estimator of (j). First, we derive the parameters 9 of 
the threshold b{t), as described in Section 4.1. Then, we use maximum likelihood 
estimator (MLE) as follows. Since the observations are independent and identically 
distributed, the log-likelihood function is given by 

n 

= ^\og fT-^irpcf)), 

i=l 

where /t~ is the pdf given by (10) with 9 replaced by 9. Then, the log-likelihood 
function can be maximized numerically to obtain the unknown parameter (j). Since 
the parameter values of p and a need to be positive, when minimizing by 

means of the function optim, we penalize negative values of p and cr by returning 
10^°. We denote by ^mle)^) the MLE of (j) derived from the threshold b with 
parameters 9. 

4.2.2. Moment estimator. A different approach consists in equating the theoretical 
moments of T^, given by Eq. (11), with the empirical moments of Tf,. In particular, 
we numerically solve a system of two equations (given by the first two moments) in 
the two unknown parameters (p = {p, a^). We denote by i^me the moment estimator 
(ME) of <p. 


0free = arg min 
0 


Estimator 

Assumption on b{t) on [To,r,] 

Unknown parameters 

Parameter conditions 

9+ 

b+{t) > b{t) 

^1 

h e [ro,r,] 


b-{t) < b{t) 

tl,t2 

To < tl < t2 < T, 

^betw 

b-(t) < bbetw(t) < 6+(t) 

none 

none 

^free 

none 

Qli A: /d2i tl 

none 


Table 1 . Proposed estimators of the parameters of the piecewise 
linear thresholds &+,&_, 6betw and bi^ee given in Section 4.1 under 
different assumptions. 
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When e is small, approximated mean and variance of Tf, are available [19, 38]. In 
particular, for a general parameter bo > Xq, we have 




vl^b) 


bo I e 
-1— exp 


'bo [ji - + 2Acr2^ 


(13) 


bo(T^ 

/i3 


2e 


+ -.{bo-l) 


fibo 


I + 2Acr2 2/1 


^0 exp 


(/i- 


2AtT2 


.(14) 


We denote by f^e moment estimator of 4> obtained from (13) and (14) when e 
is small. 


5. Simulation study. 


5.1. Monte Carlo simulations. We simulate FPTs of the Wiener process X to 
b{t) as described in [19, 28]. Applying the Euler-Maruyama scheme to the stochastic 
differential equation (1), we generate realizations of X, denoted by Xi := Ai(si), at 
discrete times Si = iAs, i>l. We set Xq = xq = 0 and As = 0.001 as time step. To 
avoid the risk of not detecting a crossing of the boundary due to the discretization 
of the sample path, at each iteration step we compute the probability that the 
bridge process s G [si,Si+i]|, originated in xt < b{si) at 

time Si and constrained to be in cci+i < 6(si+i) at time Si+i, crosses the threshold 
in between Si and s^+i. For a Wiener process, this probability is given by [15] 


P(xi,Xi+i) = exp 


2[5(si+i)^ - b{si+i){xi + Xj+i) + x,Xi+i] \ 

cr^As / 


A FPT is observed if Xi hits/exceeds the threshold b at time Si, i.e. Xi > b{si), 
or if the probability of having crossed the threshold in (si,Si+i) is larger than 
a randomly generated uniform number Ui in (0,1), i.e. P(xi,Xi+i) > Ui. In 
this case, the mid-point (s^ -I- Si+i)/2 is chosen as simulated FPT. Samples of 
size 100 are simulated for different values of (t^,A and e when 6o = 1 and /i = 
1. In particular, we consider = 0.2,0.4,1; e = 0.05,0.1,0.2,1,5,10 and A = 
0.02,0.04,0.08,0.15,0.30,0.60,1.00,3.00,5.00,10.00. These parameter values are 
chosen to cover and extend the cases of small values of e considered in [19, 38]. 
Finally, for each value of e and A, we repeat simulation of data set 1000 times, 
obtaining 1000 statistically indistinguishable and independent trials. 


5.2. Set up. In the simulations we are mainly concerned to illustrate the perfor¬ 
mance of our method under the assumption that the threshold b{t) is known, i.e. 
bo, the rate A and the amplitude e are given. Two scenarios are considered: both 
/t and are known; no information about the parameter of the Wiener process is 
given. In the first case it is of interest to evaluate the error in the estimation of 
mean, variance, CV and cdf of by comparing theoretical (11) and empirical firing 
statistics. When e is small, a further comparison with (13) and (14) is carried out. 
To measure the error in the estimation of Ft;,, we consider the relative integrate 
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absolute error (i?iAE), defined as 


RiAEiFn) = 


ihAt) - Fnmt 
E[n] 


(15) 


We replace the unknown quantities Ft^ and E[Th] by their empirical estimators, de¬ 
fined by Fn{t) = ^ X)r=i and t = X^fci Fbt/Ti, respectively. Both empirical 

quantities are based on n = 1000000 simulated FPTs, ensuring the closeness to the 
theoretical counterparts by the law of large numbers. This first scenario is meant 
to understand the goodness of our approximation through simulations. 

Another relevant question is the performance of the MLEs and MEs of n and , 
as described in Section 4.2. To compare different estimators, we use the relative 
mean error i?ME to evaluate the bias and the relative mean square error i?MSE, 
which incorporates both the variance and the bias. They are defined as the average 
over the 1000 repetitions of the quantities 


^relvM) — ; -^rel sq\f^) — 2 

/i fi 

and likewise for a^. 


5.3. Theoretical results for cdf and firing statistics of Tf,. In Fig. 3 are 

reported theoretical and empirical means, variances and CVs of Tf, as a function of 
the rate A, for small values of the amplitude e and for = 0.2,0.4 and 1. The 
given theoretical quantities are obtained from (11) for the piecewise linear threshold 
&frco- Note how the mean of does not depend on cr^, as it also happens for a 
linear threshold, while both variance and CV increase with growing We refer 
to [19] for a detailed discussion on other qualitative features of the firing statistics, 
e.g. monotonic decrease on the mean with growing A, existence of a minimum value 
for the variance, etc. What is relevant to emphasize is the excellent fit of the firing 
statistics provided by our method for any A, and for both small (cf. Fig. 3) and 
large (cf. Fig. 4) values of e. When e is small, our theoretical firing statistics 

are at least as good as those in [19, 38], outperforming them when e grows. The 

firing statistics of Tsbetw almost identical to those of while those of 

and Tb_ are slightly different for increasing e. This can be seen in Fig. 5, left 
panel, where the percentages of the i?iAE(F'T) for the four proposed estimators are 
given. As expected, the best approximation of Ft^ is provided by Ft^^ , since 
bfree is the only threshold whose parameters are obtained from a non-constrained 
optimization problem. The performance of the estimators gets worse for large 
and e. The highest error is observed for the value of A that minimizes the variance 
of Tb- However, all errors are smaller than 2%, confirming the good performance of 
the proposed estimators. 

5.4. Parameter estimation of (^, cr^). We have seen that yields the best 

approximation of Tb in terms of both cdf and firing statistics. For this reason, we 

limit our study to the estimators (p based on &free- In Fig. 6 the i?ME and the i?MSE 
of fi and (7^ are reported. As expected, the MLE provides the best estimate of p, 
while both MEs are acceptable only for small values of e. The performance of jX is 
highly satisfactory, with i?ME(A) smaller than 0.5%, and i?MSE(A) < 0.2%. Larger 
but still good i?ME and i?MSE are observed for The performance of ^mle gets 
worse for growing cr^, as shown in Fig. 7. However, except the i?ME(<5'^) for large 
values of e, all errors are between 0 and 2 — 3%. Two last remarks should be done: 





12 


MASSIMILIANO TAMBORRINO 


E(Tb) 



" o.b 2 ' o.bs' o!3 ' 1 S 6 
X 


E(Tb) 



X 


Var(Tb) 


CV(Tb) 






Var(Tb) CV(Tb) 




Figure 3. Mean (left panels), variance (central panels) and CV 
(right panels) of the FPT Tb as a function of the decay rate A 
of the threshold for small values of the amplitude e when p. = 1. 
Top panels: = 0.2. Central panels: cr^ = 0.4. Bottom panels: 

cr^ = 1. Empirical quantities from simulations (symbols), theo¬ 
retical quantities given by (11) for the piecewise linear threshold 
bfree (solid lines), and theoretical quantities (13) and (14) when e 
is small (solid gray lines). Also shown are the firing statistics of Tt 
when e = 0 (horizontal dashed lines). 


first, the i?MSE of fl for small values of e approaches the corresponding values of a^. 
Second, i?MSE(<5’^) seems not to depend on A, e and cr^, but to be equal to 2%. This 
error decreases when increasing the sample size. For example, the i?MSE(o'^) ~ 1% 
when n = 200 (results not shown). 

6. Discussion. As a consequence of the recent increasing interest towards adapting- 
threshold models for the description of the neuronal spiking activity, a need of 
suitable mathematical tools to deal with hitting times of diffusion processes to 
































HITTING TIME TO EXPONENTIALLY DECAYING THRESHOLDS 


13 


E(Tb) Var(Tb) 




CV(Tb) 



Figure 4. Mean (left panel), variance (central panel) and CV 
(right panel) of the FPT Tf, as a function of the decay rate A of 
the threshold for large values of the amplitude e when /r = 1 and 
cr^ = 0.2. Empirical quantities from simulations (symbols), theo¬ 
retical quantities given by (11) for the piecewise linear threshold 
bfrcc (solid lines), and theoretical quantities (13) and (14) when e 
is small (solid gray lines). Also shown are the firing statistics of Tt 
when e = 0 (horizontal dashed lines). 
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Figure 5. RiAEiFr^) (in percentage) given by (15) for different 
values of A and e when fi = 1. Left panel: Riae^Ft^) from thresh¬ 
old 6free (cucles), b- (triangles), &+ (rhombuses) and 6betw (gray 
circles) when e = 1 and cr^ = 0.2. Right panel: RiAEiFr^^ ) for 
cr^ = 0.2 (circles), cr^ = 0.4 (triangles) and cr^ = 1 (gray circles). 

The values of e between consecutive vertical dotted lines are fixed 
and equal to 0.05,0.1,0.2,1, 5,10, while A varies between 0.02 and 
10 . 

time-varying thresholds arises. The mathematical literature on the FPT problem is 
rich and extensive. Unfortunately, analytical solutions are not available even for a 
problem as simple (compared to others) as the one considered here, i.e. Wiener pro¬ 
cess to an exponentially decaying threshold. The closest result in this direction is 
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Rme(A) in percentage 



Rme(^^) in percentage 



X 


Rmse(p) in percentage 



Rmse(^^) in percentage 



X 


Figure 6. Dependence of RmeU^), RmseO^), Rme{^^) and 
-Rmse(o'^) (average over 1000 simulations) on A and e when X is a 
Wiener process with ^ = 1 and cr^ = 0.2. Different estimators of 
(j) = (/r, cr^) are considered: maximum likelihood estimator (/)mle 
( solid lines with triangles), moment estimator (/)me (dashed lines 
with circles) and moment estimator from (13) and (14) when e is 
small, 4>%je (gi'ay solid lines with gray circles). The values of e 
between consecutive vertical dotted lines are fixed and equal to 
0.05,0.1,0.2,1, 5,10, while A varies between 0.02 and 10. 


represented by the work of Wang and Pdtzelberger, who provide an explicit expres¬ 
sion which should then be evaluated through Monte-Carlo simulations. The idea 
behind the works of Lindner and Longtin and of Urdapilleta, was to simplify some 
mathematical difficult equations arising from the study of the FPT by linearizing 
them in e, the amplitude of the decaying threshold. As a consequence, the quality 
of the approximation rapidly decreases when e increases. 
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Rme(A) in percentage 



Rme(^^) in percentage 



Rmse(p) in percentage 



Rmse(^^) in percentage 



X 


Figure 7. Dependence of RmeO^), RmseO^), RuEif^"^) and 
-Rmse(o'^) (average over 1000 simulations) on A,e and cr^ when 
X is a Wiener process with ^ = 1 and cr^ equal to 0.2 (solid 
lines with circles), 0.4 (dashed lines with triangles) and 1 (gray 
solid lines with gray circles). Here only the maximum likelihood 
estimator 4>mle of (j) = (/i, a'^) is considered. The values of e 
between consecutive vertical dotted lines are fixed and equal to 
0.05,0.1,0.2,1, 5,10, while A varies between 0.02 and 10. 


The method proposed here has no restriction on the parameter of the thresholds 
and it is based on the simple idea of replacing the boundary by a continuous two- 
piecewise linear threshold. This allows us to derive the analytical expression of 
the FPT density to the two-piecewise threshold, and to use it to approximate the 
desired distribution. To some extent, the presence of two linear thresholds can be 
considered as a second order approximation of the problem. 

Numerical simulations show a good performance of the proposed method both 
when computing the main firing statistics, such as means, variances and CVs, and 
when calculating the FPT distribution. Different approximating thresholds have 
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been proposed. We suggest choosing the one minimizing the distance with the 
curvilinear threshold and to restrict the interval where to perform the optimization 
as described in the paper. Among the estimators of the drift and diffusion coeffi¬ 
cients of the Wiener process, we suggest applying MLE which always estimates the 
parameters reasonably well. 

The method proposed here may yield several interesting developments. First 
of all, it can be used to characterize the firing statistics of the Wiener process to 
the exponential decaying threshold, extending the previous considerations obtained 
for small values of e. Then, it may be extended to the case of a Wiener process 
with an adapting decaying threshold, as suggested in the introduction. Finally, 
our results may also be applied to all those processes which can be expressed as a 
piecewise monotone functional of a standard Brownian motion [40], as well as to 
Wiener processes with time-varying drift [19, 21]. 
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